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ABSTRACT 


In  this  paper  some  statistical  models  in  connection  with 
seismic  magnitude  are  presented.  Two  main  situations  are 
treated.  The  first  deals  with  the  estimation  of  magnitude 
for  an  event  using  a fixed  network  of  stations  and  taking 
into  account  the  detection  and  bias  properties  of  the  in- 
dividual stations.  The  second  treats  the  problem  of  estimating 
seismicity  and  detection  and  bias  properties  of  individual 
stations.  The  models  are  applied  to  analyze  the  magnitude 
bias  effects  for  an  earthquake  aftershock  sequence  from 
Japan,  as  recorded  by  a hypothetical  network  of  15  stations. 

It  is  found  that  network  magnitudes  computed  by  the  con- 
ventional averaging  technique  are  considerably  biased,  and 
that  a maximum  likelihood  approach  using  instantaneous 
noise  level  estimates  for  non-detecting  stations  gives  the 
most  consistent  magnitude  estimates.  Finally,  the  models 
are  applied  to  evaluate  the  detection  characteristics 
and  associated  seismicity  as  recorded  by  three  VELA  arrays 
(UBO , TFO , WMO) . 
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1 . INTRODUCTION 

One  of  the  main  problwuo  in  estimating  magnitude  of  seismic 

< 

events  from  network  data  is  that  the  events  are  not  always 
detected  by  all  stations.  The  standard  procedure  of  estimating 
magnitude  by  averaging  observed  magnitude  of  recording  stations 
gives  estimates  that  are  biased  upwards.  Methods  to  cope  with 
this  problem  have  been  given  by  Herrin  and  Tucker  (1972)  who 
computed  the  expected  error  introduced  by  the  averaging  procedure 
and  Ringdal  (1976)  who  developed  a maximum  likelihood  procedure. 

The  maximum  likelihood  method  given  in  this  paper  differs 

\ 

from  Ringdal' s approach  in  that  it  takes  into  account  the 
probability  that  the  event  is  detected  by  the  network.  This 

J 

will  have  effect  on  the  estimated  magnitude  for  small  events 
(events  that  are  detected  by  a small  number  of  stations).  The 
maximum  likelihood  estimator  for  magnitude  requires  knowledge 
of  the  station  detection  parameters  (threshold  and  slope  of 
detection  curve)  and  region  - station  bias. 

In  the  second  part  of  the  paper,  methods  to  estimate  these 
parameters  are  develpoed.  This  is  done  under  the  assumption  that 
the  magnitude  distribution  in  the  source  region  can  be  approximated 
with  the  usual  linear  relationship  between  magnitude  and  logarithmic 
frequency. 

The  first  method  uses  single  station  data  for  simultaneous 
estimation  of  the  slope  of  the  seismisity  curve  and  the  detection 
parameters  of  the  station.  This  maximum  likelihood  estimator  was, 
although  not  explicitly  stated,  given  by  Kelly  and  Lacoss  (1969). 
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One  interesting  result  is  that  the  distribution  of  log  (A/T) 
rather  than  magnitude  at  a station  is  independent  of  the  focal 
locations  of  the  events  and  scattering.  The  second  method  is 
developed  for  estimation  of  region  - station  bias  and  scattering 
standard  deviation.  Because  of  the  complexity  of  the  joint 
distribution  in  the  general  case  detailed  calculations  is 
carried  out  only  for  two  stations.  By  combining  analysis  made 
for  different  combinations  of  station  pairs  it  is  possible  to 
estimate  relative  bias  and  scattering  standard  deviations  for 
all  stations  in  a network. 
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A 


2 . THEORY 


2.1 


We  assume  that  the  amplitude  of  the  signal  generated  by  a seis- 
mic event  in  a specific  region  and  arriving  at  a station  can  be 
modelled  as 

(1)  A/T  = eme^eBec 
where 

A/T  is  amplitude  over  period 

m is  the  "true"  magnitude  of  the  event 

Q is  the  distance-depth  correction 

B is  the  station-region  bias,  i.e.  the  average  scattering 

effect  due  to  inhomogeneities  in  the  earth 

£ is  the  difference  between  total  scattering  and  average 
scattering . 

Taking  logarithms  of  both  sides  in  (1)  we  get 

(2)  log  (A/T)  = m + Q+  B + e 

or  letting  y = log  (A/T)  denote  the  station  magnitude  we  get 

(3)  y = m + B + Q+  e 

If  we  regard  e as  a sum  of  effects  from  many  travel  path  in- 
homogenities  or  scattering  sources  it  can  be  shown  (under  general 
assumption)  that  the  distribution  of  e is  Gaussian.  The  condi- 
tional distribution  of  y for  given  "true"  magnitude  m is  then 


(4) 


f (v/m) 


— - — expt-  ~r  (y-  (B+Q+m)  ) 7 ] 
/2tt  o 
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Here  Q is  assumed  to  be  a known  constant,  and  o2  is  the 
variance  of  the  scattering,  i.e.,  of  e.  We  will  further 
assume  that  the  station  has  a detection  curve  giving  the 
probability  that  the  event  is  seen  given  y.  We  will  here  assume 
that  this  curve  is  of  the  form 

/ \ Y 

(5)  G = / — - — exp[-  (z-G)  2 /2y 2 ] dz 

\ ' / -oo  /2tt  y 

where  G can  be  interpreted  as  the  average  threshold  and  y as 
the  standard  deviation  of  the  threshold. 


2 . 2 Distribution  of  Observed  Log  (A/T)  at  a Network  of  H 
Stations 

Let  Ni  denote  the  subset  of  station  log  (A/T)  for  given 
true  magnitude  that  are  not  seen  at  the  i:th  station.  Then,  de- 
fine the  variable  a.  such  that 

i 


a . = y . 

1 W 

if 

the 

event 

is 

seen  at  station  i; 

a . € N . 
l l 

if 

the 

event 

is 

not  seen  at  station  i 

Unless  otherwise  stated,  the  distance  - depth  corrections  (Q^) 
are  regarded  as  known  constants.  The  conditional  distribution 
of  a^/m  is  given  by 


f.Ca./m)^  -1)  da.  for  a.  * N. 

V nr 


(6) 


h.  (a  /m)  da.  = 
i i l 


/-  (Bj  + + m - G^)  ^ 


/2  , 2 

/ai  + Yi 


i for  a.  € N. 
) i i 
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If  the  scattering  effects,  i j,  1 are  i-nc*ePGndont  t^u“  j°int 

distribution  of  (a, , , . . • , a^/m)  is 

M 

(7)  h(ar  a2f...fuM/m)da1da2...daM=1E1hi(a1/m)  da. 


If  an  event  is  declared  as  detected  by  the  network  whenever  it  is 
detected  by  at  least  one  station,  v/e  obtain  the  detection 

probability  for  the  network  as 

M /-(B.  +Q.  +m-G.h  M 

(8)  P (detect/m)  = 1 - TI  — — ) = 1 - II  P(a.  € N./m) 

i=l  V / 2 . 2 ’ i=l  1 x 


/ 2 2 

!q.  + y. 

l l 


The  distribution  of  (a^,...a  /m)  given  that  the  event  is 
detected  by  the  netv/ork  is  then 


(9) 


H(a^,  a 2*  • • • , da  -j  . • . daM 


M 

1 - II  P (a.  € N./m) 
i=l  1 1 


This  is  the  conditional  distribution  of  log  (A/T)  recorded  at 
the  network.  If  the  distribution  of  true  earthquake  magnitudes 
is  known  we  can  derive  the  unconditional  distribution  of 
log  (A/T)  recorded  at  the  network.  Let  v(m)dm  denote  the 
distribution  of  earthquake  magnitude  in  a region.  The  uncondi- 
tional distribution  of  a^,...,aM  is 

CO 

(10)  g(a^ , a2» ja^daj  = j hfa^, — ,afym)da^  — daj^  v(m)dm 

—CO 

The  probability  that  all  a.  I N.  , i.e.  the  probability  of  not 
recording  is 
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‘jiving  i lif  unnond  i l i ona  1 distribution  of  I oq  amplitude  recorded 
at  the  network 


(12) 


9 (a^  # • • • i Sjuj)  da^  "** 

^1  * * * = 1-P^  € 6 Nm) 


It  is  often  assumed  that  there  is  a linear  relation  between  mag- 
nitude and  logarithmic  frequency  °f  earthquake  occurrence . 7'h.is 
corresponds  to  assuming  that  the  distribution  of  earthquake  magnitud 
is 


m; 


v(m)dm  = B exp[  G (in-m^l  ]dm  for  m > mQ 
v(m)dm  = 0 for  m < mn . 


Inserting  this  in  (12)  we  get 


j h(a1f...,a /m)da.  ..daM  exp(-Sn.)d 

J 1 n 


(14)  G(a^, . . . ,aM)da^  ...  dc^  = 


m. 


0 


M 

()  - II  P(a.  € N./m)  exp(-ftp)  dm 
i=l  1 1 


"b 


The  integrals  in  (14)  exist  for  all  m 


0 


m0  - -°°  * 


and  are  convergent  when 


In  many  cases  it  may  be  desirable  to  consider  the  distribu- 
tion of  station  magnitudes  instead  of  station  log  (A/T)  . Lett- 

ing no  denote  the  observed  magnitude  at  the  i:th  station  we  have 

(15)  m . = a . - Q . 

l l i 

and  the  distribution  of  m./m  is  then 

l 


(16) 

h *(nt  / m ) dm^h^  (m^-Kb/m)  dm^ 


m.  + Q.  - G. 

l 


* / l 

W1")  Y 


j dm.  for  m.  t N. 
/ l l r i 


r (B.  + Q*  + in  - G k 

<p  ( — i — - — — — ] for  m . £ N . 

v ' /2  2 ' 1 1 

V_  /o . + Y • 

i l 


* 

with  f^fno/m)  = f (no  + Q^^/m)  . 
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From  (lfi)  we  see  that  the  distribution  of  m./m  depends  on  the 
distance-depth  correction.  But  this  correction  is  introduced  in 
order  to  obtain  measures  of  event  magnitudes  that  are  independent 
of  the  focal  locations.  In  order  to  avoid  this  dependence  we 
can  regard  the  distance-depth  correction  as  a random  variable 
normally  distributed  with  expectation  Q_  and  standard-deviation 
Oq  . We  can  then  obtain  the  probabilities  of  detecting  an  event 
for  given  station  magnitude.  This  curve  takes  the  same  form  as 
(5).  We  get  the  detection  curve 


/m.  - G. 
l\ 


* _ * o ? 1/  2 

with  G.  = G.  - Q.  , -y  . = (YJ  + Oq  ) 


and  the  distribution  of 


(18)  hr  (nt/m)  = 


m./m  as 

i 


f . (m./m) 
l l 


* 


for 


m.  £ N. 
l l 


for  m.  € N. 
l l 


Using  h**  instead  of  h^  in  19)  and  (14)  then  gives  the  con- 
ditional and  unconditional  distribution  of  observed  station 
magnitudes . 

Comments 

We  have  here  defined  detection  of  an  event  as  seen  by  at 
least  one  station.  Using  other  possible  definitions  like  seen 
by  at  least  two  stations  will  only  chance  the  denominators  in 
(9)  and  (14),  which  in  turn  are  derived  from  Eq.  (8).  In  the 
following  sections  most  results  are  derived  using  log  (A/T) 
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as  a basis.  The  corresponding  results  in  terms  of  magnitude  are 
obtained  by  substituting  magnitude  for  log  (A/T)  and  setting 
all  distance-depth  corrections  equal  to  zero.  The  station 
parameters  G and  y are  then  interpreted  as  detection  parameters 
in  terms  of  magnitude  and  can  be  transformed  back  to  log  (A/T) 
basis  in  a way  similar  to  that  in  Eq.  (17) . 
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3.  ESTIMATION 

The  two  distrib'-LlwiiS  (9)  and  (14)  are  in  general  suited  for 
two  different  estimation  problems.  The  conditional  case  is 
useful  mainly  for  estimation  of  magnitude  of  individual  events 
while  estimation  of  structural  parameters  such  as  seismicity, 
station  bias,  etc.,  is  most  conveniently  done  in  the  framework 
of  the  unconditional  approach.  In  certain  cases,  the  conditional 
approach  can  be  used  to  estimate  structural  parameters  and  for 
those  cases  it  has  the  advantage  of  not  requiring  any  knowledge 
of  prior  distributions. 

We  shall  here  address  both  approaches  and  begin  with  the 
conditional . 


3.1  Conditional  Maximum  Likelihood 


Consider  first  the  case  of  one  station  and  one  event.  In 
this  case  the  distribution  of  observed  log  (A/T)=a  for  given  true 
magnitude  m is 


B+Q+m-G', 


(19)  H (a/m) da  = f (a/m)  $ 
or  explicit 

(20)  H(a/m)da  = — exp [ - (a- (B+Q+m ))z/2az] 

/2tt  a 


y / ^ 

— — / exp  (-t2/2)  dt / exp(-t2/2)  dt 

/2tt  — oo  / /2tF  — OO 


with  y = (a-G)/y 


x = (B+Q+m-G) //a 
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F requency 


Fig.  1 Frequency  distribution  H(a/m)  for  m = 3.8,  G = 0.0,  B+Q  = -4.0, 

o = 0.3,  y = 0.2  and  the  corresponding  Gaussian  distribution  with 
expectation  -0.2  and  standard  deviation  0.3. 

The  shape  of  the  distribution  is  shown  in  Fio.  1 together  with 
the  corresponding  Gaussian  distribution.  The  Gaussian  distribution 
corresponds  to  the  case  where  the  station  has  a probability  of 
not  detecting  an  event.  We  see  from  the  figure  that  the  effect  of 
nondetection  is  substantial.  It  is  shown  in  Appendix  1 that 


(21)  E (a/m)  = B + Q + m + - Z (x) 

with  x as  before  and 


Z(x) 


y 2n 


exp(-x2/2) 


x 

/ exp (-t2/2)  dt 
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Fig.  2 Expected  log  (A/T)  at  a station  for  given  true  magnitude.  Station 
parameters  used  are  G = 0.0,  B+Q  = -4.0,  CT  = 0.3,  y = 0.2. 


Figure  2 shows  E(a/m)  for  B+Q  = -0.4,  G = .0,  a = 0.3,  y = 0.2. 

As  can  be  seen  from  the  figure,  the  bias  is  large  for  events  of 
smaller  magnitude  and  also  quite  large  for  events  around  the 
detection  threshold. 

In  this  case  it  is  readily  shown  that  the  maximum  likelihood 
estimate  corresponds  to  solving  the  following  equation  for  n 

1221  a = B + Q + m + - 2 Z(x) 

/^TT2 


— 
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In  the  case  of  a network  consisting  of  several  station, 
one  might  estim.-.  t®  the  magnitude  for  each  station  that  detected 
the  event  according  to  the  above  approach  and  then  average 
these  values  to  obtain  an  estimate  of  the  true  magnitude. 

However,  if  the  stations  that  did  not  see  the  event  were 
operating  properly,  the  information  that  some  of  the  stations 
did  not  see  the  event  is  useful  information  and  should  be 
utilized  in  the  estimating  procedure.  This  was  first  done 
by  Ringdal  (1976)  who  considered  the  case  with  y = 0,  i.e., 
the  stations  detect  the  event  as  soon  as  the  'observed 
magnitude'  nu  is  larger  than  G.  His  results  will  differ 
from  those  given  here  as  he  included  in  the  sample  space 
the  cases  where  the  event  was  not  seen  at  any  of  the  stations. 
This  will  have  effect  on  the  estimated  magnitude  for  smaller 
events  but  for  larger  events  and/or  large  number  of  stations 
the  differences  will  be  small.  The  reason  for  this  is  that 
the  probability  of  not  seeing  an  event  tends  to  zero  as  the 
true  magnitude  increases  and/or  the  number  of  stations  increases. 

Suppose  that  we  have  a network  of  M stations  recording  an 
event.  Suppose  further  that  the  event  is  seen  by  at  least  one 
station.  The  corresponding  log  likelihood  is 

M M (B.  + Q.  + m - G.) 

(23)  log  L = Z log  h.  (a./m)  - log(l  - n * ( — j-— _ 1 

i=l  11  i=l  V /2  2 

/oi  + Yi 
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We  can  see  ucu;  that  the  likelihood  is  not  a product 
of  'individual  likelihoods'.  If  we  had  considered  the  case 
with  events  seen  at  all  stations  or  a specific  station, 
the  likelihood  (19)  would  have  been  a product  of  'individual 
likelihoods ' . 

It  is  shown  in  Appendix  1 that  the  likelihood  estimator  for 
magnitude  is  consistent  and  its  asymptotic  standard  error  is  otiven. 

Can  we  use  the  conditional  approach  to  estimate  the  parameters 
B^,  , ck  and  y^  and  magnitude  simultaneously?  To  try  to  an- 

swer this  question  we  assume  that  we  have  a sample  of  N events 
such  that  each  event  is  seen  and  recorded  by  at  least  one  station. 
In  terms  of  "true'-  log  (A/T)-Q  arriving  at  the  stations  we  have 
the  following  model 


(24) 


il 


B . + m . + e . . 
i 1 il 


i=l , 2 , . . . ,M;  j=l , 2 


/ N 


If  we  had  had  the  case  where  every  event  was  alwavs  seen  at  every 
station  this  would  have  been  the  standard  model  for  analysis  of 
variance  of  a two  way  classification.  In  order  to  make  the  model 
identified  we  have  to  impose  a normalizing  condition  on  the  B^rs 
(or  m^ts).  The  reason  for  this  is  that  these  unknowns  are  iden- 
tified only  up  to  an  additive  constant.  We  will  here  adopt  the 
M 

condition  T.  B.=0,  that  is,  the  average  station  bias  is  set 
i=l  1 

equal  to  zero. 


In  the  analysis  of  variance  model  with  fixed  effects  the 
likelihood  does  not  have  a maximum  if  we  treat  the  o.:s  as  un- 

l 

known.  This  will  clearly  also  be  the  case  when  we  have  missing 
observations,  i.e.  some  stations  fail  to  renort  an  event  bacause  of 
the  detection  properties  of  these  stations.  However,  even  if  we 
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treated  the  o^.  :s  a s known  constants  there  are  difficulties  be- 
cause the  nun.’.  ~r  of  unknowns  increases  with  sample  size.  That  is, 
as  the  number  of  events  increases  the  number  of  unknown  event 
magnitudes  increases.  Therefore,  the  standard  methods  for  inves- 
tigating the  properties  of  the  maximum  likelihood  estimator  do  not 
apply.  The  reason  for  this  is  that  the  estimator  for  itk  is  not 
sufficient  independently  of  B^ , and  y ^ • And  as  we  cannot 


obtain  a simple  relation  expressing  the  estimator  for 


m . 
1 


in 


terms  of  data  and  the  parameters  B^,  and  y^  we  cannot  use 


the  approach  used  by  Christoff ersson  (1970)  in ^connection  with 


estimating  factor  loading  in  the  case  of  missing  observations. 


So,  at  least  for  the  time  being,  we  have  to  regard  the  condi- 
tional approach  as  a heuristic  one  for  estimating  the  station 
parameters . 


For  the  case  of  just  estimating  the  magnitude,  the 
likelihood  estimator  seems  to  be  fairly  insensitive  to  moderate 
changes  in  the  present  parameters.  However,  it  is  important 
that  the  stations  that  do  not  report  an  event  have  been  operating 
properly  at  that  time.  Otherwise,  there  will  be  large  biases 
in  the  estimated  magnitudes.  The  asymptotic  standard  errors, 
on  the  other  hand,  are  directly  related  to  the  present  level 
of  scattering,  i.e.,  on  the  values  and  also  to  some  extent 
on  the  y^:s.  Tables  1 and  2 and  Fig.  3 show  an  application  for 
a simulated  network  consisting  of  15  stations  to  an  aftershock 
from  Japan.  Here  all  are  preset  to  0.3  and  all  y^  to  0.2. 

In  the  figure  there  are  some  outliers.  The  first,  event  no.  2, 


16 


Fig.  3 Average  and  maximum  likelihood  m^  for  a simulated  network 
of  15  stations.  Aftershock  consisting  of  72  events. 

illustrates  the  importance  of  the  earlier  mentioned  condition 
that  the  stations  have  to  be  in  operation  and/or  that  no 
extremely  high  noise  levels  or  interfering  events  are  present 
because  this  event  is  so  large  that  it  should  have  been  seen 
by  all  the  stations.  Therefore  it  is  probably  best  to  have 
a variable  threshold,  for  example  related  to  the  noise  level 
in  the  time  interval  preceding  the  signal,  as  suggested  by 
Ringdal  (1976).  The  other  outliers,  events  3,  4,  40  and  47, 
are  small  events  seen  by  3,  1,  1 and  1 stations,  respectively, 
and  reflect  the  bias  in  the  conventional  method  for  magnitude 
estimation . 
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TABLE  1 

ESTIMATES  OF  STATION  BIAS  Bi#  AND  STATION  DETECTION 

THRESHOLD  G. 

1 


Station 

B. 

1 

Qi 

Gi 

LAO 

0.07 

3.7 

-0.1 

MBC 

0.29 

3.9 

+ 0.2 

NAO 

0.00 

3.7 

0.0 

RES 

0.38 

3.7 

+0.5 

HFS 

0.00 

3.7 

+0.1 

UBO 

-0.16 

3.7 

+ 0.1 

KBL 

0.09 

3.9 

+0.2 

FFC 

-0.03 

3.7 

+0.6 

BLC 

0.20 

3.7 

+ 0.9 

ALE 

0.05 

3.7 

+0.8 

FBC 

0.05 

3.7 

+ 0.8 

CHG 

-0.39 

3.7 

+ 0.3 

COL 

-0.21 

3.8 

+0.4 

FCC 

-0.15 

3.7 

+0.7 

YKC 

-0.21 

3.7 

+ 0.7 

STATION 
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3.2 


i 

Unconditional  Maximum  Likelihood  Estimation 

Consider  first  t!  " _se  with  just  one  station.  From  « 

« 

Eq.  (14)  we  get,  letting  -*■  -°°  j 

CO  / oo  I 

(25)  G(a)da  = / h(a/m)da  exp(-Bm)  dm y j [l-P(aGN/m)]  exp(-Bm)  dm  ‘ 

Now  J 


(26)  / h(a/m)  exp(-gm)  dm  = / 


-°°/2n  a 


exp [ - (a- (B+O+m) ) 2 /2a7  ] 


exp(-Bm)  4’^-^^jdn  = const.  exp(-Ba)  4) 

< 

OO  OO  ' V 

(27)  / [1-P(a  N/m) ] exp (-Bm)  dm  = / <j>  | B_+Q+m — G j eXp(_gm)  = 

-oo  -oo  \/a2  +y 2 / 

= const. 


Thus,  the  distribution  of  seen  events  at  the  station  is 
proportional  to 

(28)  exp(-Ba) 

Evaluating  the  proportionality  constant  (see  Appendix  2)  we 
find  that  the  distribution  is 

(29)  G(a)da  = exp (BG-y2 82/2)  exp(-Ba)  $ 

We  note  that  this  distribution  is  independent  of  the 
scattering  standard  deviation,  the  region-station  bias  and 
distance-depth  correction,  i.e.,  the  observed  distribution 
depends  only  on  the  seismicity  parameter  8 and  the  station 
detection  parameters  G and  y • 


ci""G 


da 
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Turning  to  the  likelihood  estimator  we  have  for  a sample  of  N 
independent  observaLiuns,  a^,  a2» . • • • > 

12  2 N N /ai  ~ G\ 

(30)  log  L = N log  B + N S G - i N -3  l a.  + E log  $ I— , 

* i=l  i=l  v / 

The  derivatives  of  log  L with  respect  to  3/  G,  and  Y are 


(31) 


(32) 


(33) 


6 log  L 

6 3 


N 

= £ + N G - N 

p 


N 

I 

i=l 


a. 

x 


6 loq  L ...  1 ” *’  (yi) 

h ' 6 y 


N 


<5>'  (yi) 


6-  = N 32  Y - i Yi  * (Yi) 


where  y. 

J l 


The  likelihood  estimator  is  defined  as  the  solution  to 


6 log  L _ 6 log  L _ 6 log  L 
6 3 6 G 6 y 


which  corresponds  to  the  maximum  of  log  L. 

Although  not  explicitly  stated  these  likelihood  equations 
(eq.  31-33)  were  obtained  by  Kelly  and  Lacoss  (1969)  who  in  addi- 
tion were  estimating  the  total  number  of  earthquakes  in  a given 
time  period.  The  likelihood  estimator  for  3 considered  by  Aki 
( 1965)  can  be  obtained  as  a special  case  from  the  above  equations 
by  putting  y=0  and  G=min (ax , . . . ,aN) . 
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To  illustrate  the  above  method,  data  in  the  distance 
range  30°-60°  observe;'  the  stations  UBO,  TFO  and  WMO  have 
been  analyzed.  The  results  in  terms  of  base  10  logarithms  are 
shown  in  Figs.  4,  5 and  6.  From  the  figures  we  see  that  the 
model  fits  the  data  reasonably  well.  However,  for  UBO  (Fin.  4) 
there  is  statistically  significant  deviation  between  model  and 
data.  On  the  other  hand,  the  number  of  observations  is  large 
(4562)  so  that  the  probability  of  detecting  small  deviations 
is  large.  The  maximum  deviation  between  data  and  model  is  0.02. 
Whether  this  deviation  is  of  any  practical  significance  will 
depend  on  the  situation  in  which  the  model  is  applied.  For  example, 
in  seismic  risk  studies  this  deviation  may  be  of  very  great 
importance . 


Frequency 


Fig.  4 UBO  distance  30  -90  . 4562  events.  Observed  and  theoretical 
l°gi0 (A/T)  cumulative  distributions.  0 = 0.85  (.017), 

O = 0.19  (.010),  Y = 0.11  (.007).  Standard  errors  within 
parenthesis . 


- 2 0.0  Q2  04  06  08  10  12  14  16  18  2.0  2.2  24  2.6  2.8  3.0 


Fig. 


5 TFO  distance  30°-90°.  3725  events.  Observed  and  theoretical 
log  (A/T)  cumulative  distributions.  3 = 0.93  (.019), 

O = 0.09  (.008),  y = 0.10  (.005).  Standard  error  within 
parenthesis . 


j < ’ifUfnry 


WMO  distance  30°-90°.  2321  events.  Observed  and  theoretical 
log  (A/T)  cumulative  distributions.  3 = 1.16  (.035), 

O = 0.40  (.015),  y = 0.16  (.007).  Standard  error  within 
parenthesis . 


Fig.  6 
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Looking  at  the  estimated  8-values  we  see  a rather  larqe 
variation,  $UBQ  = 0.85,  &mQ  = 1.15  and  ^Tp0  = 0.93.  This 
difference  may  be  attributed  to  the  different  seismic  area 
(see  Figs.  7,  8 and  9)  that  are  sampled.  The  southern  part 
of  Japan  is  outside  the  distance  range  30°-90°  for  WMO  which 
may  explain  the  relatively  high  6-value  for  WMO.  On  the  other 
hand,  it  has  been  suggested  that  the  linear  seismicity  model 
is  valid  only  for  a limited  magnitude  range  and  that  the  slooe 
of  the  seismicity  curve  increases  with  magnitude.  The  relatively 
high  detection  threshold  and  8 for  WMO  is  supporting  this  latter 
hypothesis . 

Whereas  the  magnitude  distribution  for  one  station  can  be 
used  to  obtain  reliable  estimates  of  the  detection  parameters 
G and  y no  information  about  the  station  bias  B and  the  scattering 
a can  be  obtained  as  the  distribution  of  observed  magnitude 
at  a station  is  independent  of  these  parameters.  To  estimate 
these  parameters  we  would  like  to  consider  the  distribution 
(14)  for  the  general  case  of  M stations. 


r 


no  T0N70  FOREST  HK. 


Fig.  9 Azimuthal  projection  with  center  at  TFO.  The  distance  curves 
from  origin  are  30°,  45°,  60°,  75° , 90°.  Azimuths  are  0,  30, 
60,  90  a.s.o.  degrees. 
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Unfortunately,  evaluation  of  the  integrals  in  (l'1)  for 
the  general  case  lead  J very  complicated  expressions  which 
do  not  seem  to  be  ust_-lul  from  a practical  point  of  view 
(except  of  course  for  the  case  with  events  seen  jointly  by 
all  stations,  but  then  it  is  very  difficult  to  obtain 
enough  events  for  analysis) . Using  numerical  integration 
does  not  seem  to  work  with  more  than  a few  stations. 

However,  for  the  case  M=2  and  events  seen  jointly  by 
the  two  stations  the  unconditional  distribution  is  not  too 
complicated.  In  this  case  the  distribution  turns  out  to  be 
(see  Appendix  2) 


35) 


where 


G(a  a ) =====-J 

1 2 /a2+a2/27  V Yi 


ai"Gi\  /a_-G„ 


’1\  (2  2\ 

-j  exp  (-y/2) 


(*  (c^)  exp(«2^  + 

$ (a3)  exp  («4>  ) 


(36)  y = (a^a^B^-B^))2  + 3 (2o2  (a^-Q^  + 2a2 (a 2-B2~Q2 ) - 


- a2  a2  6) ) / (a2  + a2) 


(37)  rX;L  = ((B1-KJ1  - B2-Q2)  - (Gx-G2)  - (a2  + y2)  3) 


, 2 , 2 2 2 

al  + °2  1 +Y2 


(38) 


ou  = 


3(b2+Q2"G2)  + 62(o2  + Y2)  / 2 


(39)  a 3 = ((B2+02  - Bj-Q^  - (G^)  - (a2  + Y2)  3) 


A 


,2  2 2 2 
°1  + °2  + Y1  +Y2 


<4  = 6(B1  + Q1  - Gx)  + 62(a2  + y2)  / 2 


(40) 
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When  comparing  two  stations  it  is  a common  practice  to  plot  the 
magnitude  at  one  station  for  given  magnitude  at  the  other.  It 
is  shown  in  Appendix  2 that  the  conditional  distribution  of 
log  (A/T)  at  station  2 for  given  log  (A/T)  at  station  1 is 


(41)  g(a2/a1)  = 


a2"G2) 


— - — exp[  - (a2-a^-  (C^-Q^)  -B) 2/2o2 ] 
/2tt  a 


a1+B+(Q2-Q1)-G2 


/^T 


where 


(42)  B = B2  - BL  - 6 

(43)  a2  = a2  + a2 


The  parameters  that  are  identified  are  B,  G2,  a2  and  y2,  and 
these  can  thus  be  estimated.  The  corresponding  likelihood 
estimator  is  given  in  Appendix  2. 


This  distribution  is  of  the  same  type  as  the  distribution  of  ob- 
served log  (A/T)  for  given  true  magnitude  (eq.  (20))  so  the  first 
two  moments  of  (41)  are 


(44) 


E(a2/ai) 


EH-(Q  -qi)+a1+  —a- 
1 & 


+r; 


<P'  (x) 
4>  (x) 


and 


(45) 


D^(a>/a^)  = 


4 

o 


(aW2) 


*'  M („  + $l(x)\ 

4>  (x)  \X  4>  (x) ) 


U 
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where  B and  a given  by  (42),  (43)  and 


(46) 


x = 


VB+(Q2-^'  ~°2 


If  we  had  also  considered  events  not  seen  at.  station  2 (but  seen 
at  station  1)  we  would  have  obtained  the  following  distribution 


I 

(47)  g* (a2/a1)da2  = j 


— - — exp[-  (a?-a1-(01-Q„)-B)/2a2]4>.  da2 

/Ho  2 112  ^ y I 


if  seen  at  st.  2 


$ 


- (3.,+B  (Q2-Qj_)  -G2); 


\ /o 2 +y  r 


if  not  seen  at  st.  2 


B and  o are  given  by  (43) . This  distribution  is  of  the  same 
type  as  (6 ) . 

Ringdal  (1975)  based  a likelihood  estimator  for  the 
detection  parameters  on  the  probabilities  of  detecting  an 
event  at  station  2 for  given  observed  magnitude  at  station  1, 
i.e.,  essentially  on  Eq.  (47).  However  with  this  latter 
approach  it  is  not  possible  to  separate  the  effects  of 
scattering  (a2)  and  slope  of  the  detection  curve  ( y 2 ^ • 


i 


4 


1 


< 

i 
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To  illustrate  the  method  based  on  the  distribution  (41)  the  same 
data  as  for  the  sinyle  station  case  has  been  used.  The  analysis  has 
been  made  for  both  log  (A/T)  and  magnitude  data.  Letting  denote 

the  parameters  in  the  model  with  magnitude  data,  the  relation  to 
parameters  with  log  (A/T)  data  is 


(48) 


B*=  B ; 


G2+Q2J 


* 

a = 


a 


2 

VSQ2 


■N 


2 

where  Q is  the  average  distance-deepth  correction  and  is  the 

corresponding  variance. 


TABLE  3 

LOG  (A/T)  DATA 


Station 

Ref . station 

B G 

S 

Y 

UBO 

WMO 

-.04  .07 

(.01) (.01) 

.37 

(.01) 

.05 

(.01) 

TFO 

.15  .11 

(.01) (.02) 

.37 
(.01  ) 

.08 

(.02) 

UBO 

-.34  .37 

( .02) ( .02) 

.37 

(.01) 

.17 
(.01  ) 

WMO 

TFO 

-.07  .39 

(.02) (.02) 

.35 

(.01) 

.18 

(.01) 

UBO 

-.44  .03 

(.01) (.01) 

.37 

(.01) 

.08 

(.01) 

TFO 

WMO 

-.29  .01* 

(.01 ) (.02) 

.36 

(.01) 

.05 

(.02) 

* 

Lower  limit  for  G. 

Standard  error  within  brackets. 
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TABLE  4 

MAGNITUDE  DATA 


Station 

Ref .station 

B 

G S 

Y 

Q 

SQ 

WMO 

.14 

(.01) 

3.88  .36 

(.03) (.01 ) 

.24 

(.02) 

3.71 

.16 

UBO 

TFO 

-.05 
(.01  ) 

3.82  .37 

(.03) (.01 ) 

.21 

(.02) 

3.72 

.17 

UBO 

-.32 

(.02) 

4.17  .37 

(.04)  (.01) 

.24 

(.01) 

3.75 

.17 

WMO 

TFO 

-.06 

(.02) 

4.20  .35 

(.04)  (.01  ) 

.26 

(.02) 

3.74 

.16 

UBO 

-.42 

(.01) 

3.77  .36 

(.03) (.01 ) 

.22 

(.01) 

3.71 

.17 

TFO 

WMO 

-.26 

(.01) 

3.67  .35 

(.03)  (.01) 

.20 

(.03) 

3.71 

. 1 5 

Standard  errors  within  brackets. 
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The  results  of  the  analysis  are  shown  in  Tables  3 and  4.  In 
Fig.  10  a plot  of  the  observed  average  magnitude  at  '.VMO  for  given 
magnitude  at  UBO  is  shown  together  with  the  expected  average 
magnitude  calculated  from  Eg.  (44).  There  were  no  difficulties 
in  solving  the  likelihood  equations  except  for  one  case;  using 
log  (A/T)  data  in  estimating  the  parameters  for  TFO  with  WMO 
as  reference  station.  This  problem  is  probably  caused  by  in- 
sufficient coverage  of  events  in  the  neighborhood  of  the  detection 
threshold  for  TFO  due  to  the  high  detection  threshold  for 
the  reference  station.  Anyhow,  by  restricting  the  oermissible 
range  of  variation  for  the  estimated  threshold  so  that  it  had 
to  be  larger  than  the  smallest  observed  log  (A/T)  + 0.01  it 
was  possible  to  obtain  reasonable  estimates  of  the  parameters . 
Comparing  the  results  based  on  log  (A/T)  and  magnitude  we  find 
close  agreement.  If  we  use  Eq . (47)  and  transform  log  (A/T) 
parameters  to  magnitude  or  vice  versa,  the  results  will 
differ  only  slightly.  There  is  a tendency  for  magnitude  data 
to  give  somewhat  larger  numerical  values  for  both  G and  y. 
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In  the  distribution  (41)  it  is  neighter  possible  to  identify  the 
individual  scattering  variances  nor  the  region-station  biases 
But  if  we  can  assume  that  the  parameters  are  at  least  approximatly 
constant  over  the  different  regions  we  can  combine  the  results 
for  all  possible  combinations  of  station  pai>-3  and  obtain 
estimates  of  scattering  variance  and  region  station  bias.  Using 
log  (A/T)  data  we  have  for  the  different  combinations  (see  eq.  43). 


:2  + 
UBO 

_2 

WMO 

= 

.37 

2 

UBO+ 

2 

°TFO 

= 

.37 

2 

UBO+ 

2 

°WMO 

= 

.37 

2 

WMO 

2 

°TFO 

= 

.35 

.2 

UBO+ 

2 

°TFO 

= 

.37 

,2  + 
WMO 

2 

TFO 

= 

.36 

Solving  this  system  by  least  squares  we  obtain  the  following  estimates 
of  the  scattering  standard  deviations:  orT  = .27  ; aiwrt=  .25  and  o_  = .25 
Similarily , for  magnitude  data  we  get  aril  = .27  ; a.  ==  .25  and  a = .24 

UdU  WMU  lrU 

Turning  to  the  region-station  bias  we  obtain  from  eq.  (43)  using 
log  (A/T)  data 


bubo“ 

bwmo~  e 

°UBO- 

1 

O 

4* 

bubo" 

btfo"  b 

»UBO* 

.15 

bwmo“ 

bubo"  b 

2 

WMO~ 

-.34 

bwmo“ 

btfo"  6 

2 

WMO" 

i 

• 

o 

btfo“ 

bubo“  b 

2 

TFO- 

-.44 

btfo" 

bwmo"  b 

Q 

•-3  to 
'tl 

O 

II 

-.29 
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Wt  first  note  that  the  region-station  coefficients  can  only  be 
determined  up  to  an  unk.  additive  constant.  We  therefore  adopt 
the  normalizing  condition  0.  If  in  addition  we  substitude 

IHJ 

2 

the  estimated  scattering  variances  obtained  from  (49)  for 
°WMO  and  CTF0  We  Can  SOlVe  (50)  f°r  BUBO'  BWMO  and  6 by  leaSt 

A A /V 

squares.  This  gives  BUBo=  ' BWM0=  • ^ and  2.78, 

(in  base  10  logarithms  8=  1.20).  This  estimate  of  B must  be  regarded 
as  very  unreliable  compared  to  the  estimates  of  ByB0  and  BWMO 
because  the  diagonal  element  in  the  inverse  of  the  moment  matrix 
obtained  when  solving  (50)  is  more  than  100  times  the  elements 
corresponding  to  BUBQ  and  BWMQ.  Using  magnitude  data  we  obtain 
BllBO=  ' BWMO=  • 1 2 and  8=  2.53  (in  base  10  logarithms  8=  1.10). 

The  reliability  of  this  B-value  is  of  the  same  order  as  for  log  (A/T) 
data . 

It  must  be  noted  that  the  results  concerning  bias  and  scattering 
are  based  on  the  assumption  that  the  parameters  are  the  same  for  all 
sampled  regions  and  that  the  analysis  of  single  station  data  in  licated 
that  the  slopes  (8)  may  not  be  equal.  This  difficulty  can  be  overcome 
by  using  data  from  one  specific  region.  However,  in  this  paper,  data 
is  used  only  to  illustrate  the  statistical  methods.  And  for  this 
purpose  data  is  considered  to  be  of  sufficient  quality. 
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Appendix  1 


The  condit-.. distribution  of  observed  log(A/T)  at  a static 


given  by  e^.  ( ly , 


(Al.l) 


H(a  /m)  = exP[-(a-(B+G+m))2/2a2]  <J>(y) / $(x) 


where  y = (a  -G)/y;  x =(  B+Q-Hn-G)/ /o2+y? 


and  $(x)  = — / exp(-t2/2)  dt 

/2rr  _oo 

To  obtain  the  expectation  of  observed  magnitude  for 
magnitude  we  first  evaluate 


given  true 


(A1.2) 


F = f ( a -(B+Q+m))  exp[-(a-(B+Q+m)  )2/2o2]  <j>(-SlG-) 

_oo  Y 


Put  u - a - ( B+Q+m ) and  G = B+Q+m-G 

o 


(A1.3) 


1 °°  u+G 

F = /2^o  1 u exP (~u2/2ct2)  4>( du 

— oo  y 


{Al.U) 


3u  ex^  u2/2o2)  - -u/a2  exp(-u2/2a2) 


(A1.5) 


F = [ j==-  exp(-u2/2 a2)  4>( &)  ] 

VtrTT  V —a 


+ f ^ exp(-u 2 /2a2)  — i—  exp[-(u+Gn)2/2Y2]  du 

--  Y ° 


o + / — expt-(u2Y?+(u+G  )2o2)/2o2y2] 

_no  C ' ^ 


(A1.6) 


uV  * <wo)V  = ^ + Sji  2 
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We  then  find 


(A1.7)  F = exp[-G2/(cr  ■ ,'2)2]  / ~~  exp[-(a?+y2)  (u+Goa2/(a2+Y2 ) )2/2<j2y  ']  du 

00  c"  ^ 

and  finally 

1*1-8)  F - »^[-°S«("2^2)) 

Then 

(A1.9)  E(u/m)  = F/4>(x) 

and 

(A1.10)  E(a  /m)  =B+Q+r>F/4>(x)=B+Qhnt^=7 

B-KHm-G 

where  x = ^77^ 

To  obtain  the  variance  of  observed  magnitude  for  given  true  magni- 
tude we  first  evaluate 


(Al.ll) 


1 

> 


OO-I  9 U+'J 

F = / rr=~  -7  (^7  -l)  exp(-u2/2n2)  $(— — -)  du  = 

/2iva  a a Y 

.00 

u+G  00  00 

= l^oT  eXp{"u2/2o7)  t(T^)]  + / ^7  exp(-u2/2a2) 

— 00  — CO 

exo[-(u+0  ) 2/2y2  j du 
o 

= / -u— ^ exp(-u2/2a2)  exp  [-(u+G  )2/2y2]  du  = 

27raJY  o 

mmOO 


-Gn  exp[-02/2(a2+Y2 ) J 

y^iT  /cj2+y^  (o2+y2) 


(o2+y2  ) 


<t>'(x)  with  x as  before. 


This  gives 

(A1.12)  E(u2/m)  = a2  - (0^2)7(x) 


-inl- 


and 

(A1.13)  D2 (u/m)  = E(u2/m)  - E2(u/m)  = a2  { + ) 

(a^+Y  v'Kx)  4>(x ) 

f 

Thus,  the  variance  in  the  observed  distribution  is 


(A1 .lU ) D2( a /m) 


o4  (x) 
(o2+y2 ) $(x) 


(x  + 


(x) 

0(x) 


) 


It  can  be  easily  shown  that 


(A1.15) 


lim 

X-VOO 


4>’  (x) 
*(x) 


(x  + 


$’(x) 

4>(x) 


) = 0 


and 


(A1.16) 


lim  $ 1 (x ) 
x-*— 00  $(x) 


(x  + 


$ ' ' ( X ) 

4>(x) 


) = 1 


This  gives 

(A1.1T ) lim  D2(a/m)=o2 
m->°° 

and 


(A1.18) 


Turning 

(A1.19) 


iS-  D2U/m)  = o2 

to  the  likelihood  estimator  we  have 
log  L = log  H(.a  /M)  = 

= const.-(a-(B+Q+m ) )2/2o2  - log  4>(x) 


(A1.20)  3 log  L 


3m 


= ( a _ ( B+Q+m ) ) / o - 


1 $jjx) 
/o2+ Y2  <t(x) 


Solving  s L = 
0 3m 


0 then  leads  to  solving 


(A1.21) 


a 


= B+Q+m+ 


o2 

/oz+y2 


*'(x) 

4>(x) 


1»2  - 


I 


i 

i 


or  letting  z(x) 


$'  (x) 
♦ (x ) 


(A1.22)  a = B+Q+K+  ^7^7-  z(x) 

Turning  to  the  case  of  several  stations,  the  likelihood  is  from 
(eq.  19) 

M M 

(A1.23)  L(m)  = n h.  (a  . /m)/(l-tt  P(a.e  N./m)  = 

.,11  . , i i 

i=l  i=l 

= L*(m)/P*(m) 


with 


M 

(A1.2M  L*(m)  = tt  (h.(a./m) 

i=l  1 1 

M 

(A1.25)  P*(m)  = 1 - 17  P(s.£TI./m) 

i-1  1 1 


r=—  exp[-(a.-(B.+Q.+m) )2/2o?3 
f v2vo±  ill  l 


if  the  event  is 
seen 


(A1.26)  h (a  /m) 


— (m+B-j+Q1  — G 1 

V <t>(  / 2+  7 r ~ ~ ) if  the  event  is  not  seen 


and 

-(m+Bi+Gi 

(A1.27)  P(a.e»,/m  = » , g.- g * 

1 1 *ai  Yi 

It  was  mentioned  earlier  that  this  likelihood  is  not  a product  of 
independent  likelihoods,  so  the  usual  asymptotic  results  for  maximum 
likelihood  do  not  apply  directly.  However,  by  showing  that  the  maximum 
for  L(m)  is  the  same  as  for  I.*(m)  when  M tends  to  infinity  we  can 


- 1+3  - 


apply  the  asymptotic  formulaes  to  L*(m)  which  is  a regular  likelihood. 

We  have 

(A1.28)  log  L(m)  = log  L*(m)  - log  P*(m) 

Provided  ,cri'Y^  are  well  behaved  for  all  i,  that  is,  are  such 

that  each  station  has  a non  zero  probability  to  see  the  event  it  follows 
that  log  L(m)  tends  to  log  L*(m)  for  every  finite  m.  Specifically 
assume  that  there  exists  a <$>0  that 

_(ai+Bi+Qi~Gi 

(A1.29)  6 £ 4>  ( /~T  ■ 2 1-6  for  all  i 


(A1.30)  1 - (1-6  )M^.  P*  (m  ) £ 1 

o 

And  as  M tends  to  infinity  we  have 

(A1.31)  liffl  P*(m  ) = 1 
M-x»  0 

It  then  follows  directly  that 
(A1.32)  lim  log  P*(m)  = 0 

n 

(A1.33)  9 3mP~P*'^'  = ° f°r  a11  finite  m and  n=1>2»3 

From  this  it  follows  that  the  asymptotic  properties  of  the 

estimator  obtained  by  maximizing  L(m)  is  the  same  as  that  obtained 
by  maximizing  L*(m).  As  the  latter  is  a regular  maximum  likelihood 
estimator  we  have,  letting  mQ  denote  the  true  magnitude  of  the  event 
and  m the  likelihood  estimator,  that  (m-mQ)  tends  to  a normally 


- >t)( 


i i ..',r  i iK'-ii  variaM'  wi  th  •>-tati'»n  /.>ri<  and  variant.1 


(A1.3U 


) D2(4i  (m-m  ))  = [(£  E.(— 1o-^  L*(m) 


^ E(^_u_w!w)/n=Tn  )2)] 

M dm  o 


M 

(A1.35)  1 E[(91o|L*(m)  ; )2]  = 1 z E[  (11°£  hi (alM ) j m=m  )2; 


o J M 

i=l 


^ -i  ■ ( B +Q  ^ -fm„ ) 

1 j- — 1 — ° if  the  event  is  seen 

at 


(A1.3C) 


3 log  h. (a./m) 


/ai+n 


*/a?+yT 


(-Ji!Ia±Eij 

/opyT 


if  the  event  is 
not  seen 


or  letting 


(A1.37) 


(A1.38) 


we  have 


u.  = a.  - (B.+Q.-ttn  ) 
; 1 1 O 


x.  = (m  +B.-KJ.+G. 

1 Oil1  1 i 


(A1.39) 


31og  h-;  (a-j  /m) 


/ m=m 


u.  /a? 
i i 


■ t'(xj  ) 
/o|+y?  't(-Xi) 


(A1.U0)  E((-  ^ (a4/m ) i m=m  ^2)  = / -==  '■£  exp(-u2/2o?)  du.  + 


/2 no 


+ 7[+y7  (I(^J')2  ^(-xi} 

From  eq.  (Al.ll)  it  follows  that 


- U5  - 


(Al.Ul) 


And  finally 
(A1.U2) 


E((3  log  hi(Vm)  / m=m  )2) 
om  o 


= *u.)  4 + (^4-  *.)  = b. 

i a.  or+y.  $(-x. ) 1 i 


i i i 


i M 

Di 2(*ft(m-m  ))  = l/(jj  Z b . ) 
O M . , l 

1=1 


- 1*7  - 


Appendix  2 

The  unconditien- i 'istribution  of  observed  log(A/T)  at  a station 
is  from  eq.  (23)  proportional  to 

( A2 . 1 ) G*(a)  = exp (-6a)  4’(^--i) 

As  the  integral  of  the  distribution  equals  1 we  can  determine  the  proportionality 
constant  by  evaluating 

( A2 . 2 ) / exp(-Ba)  ♦(^-^)  da  = [— | exp(-Ba)  4>(^-^-)] 

— oo  Y p Y 


+ — f exp(-Ba)  — — — ex->[-(a-G)2/2y2]  da 
6 -»  /2v  y 

OO 

= 0 + x J £=“  exp[-(a-(G-2By) )2/2y2]  exp[-BG+B2y2/2]  da 

D vdTiy 

.oo 

= i exp(-BG+B2Y2/2) 

8 


Thus,  the  distribution  of  observed  log  (A/T)  is 

(A2.3)  G(a)  = B exp(  BG-B2y2/2)  exp(-Ba)  ll>(^i) 

Suppose  we  have  a sample  of  N independent  events  recorded  at  the 
station.  The  likelihood  estimator  of  the  parameters  B,G,y  are  obtained 
by  maximizing  the  logarithm  of  the  likelihood.  We  have,  letting  a^ 
denote  the  observed  log  (A/T). 


FRBOkDlNQ  PACK  BUMUWOT  FU*£i> 


U8  - 


(A2.U) 


N 

log  L(8,G,y^=  7 log  G(a^)  = 

1 N 

= N log  p + N8G  - i Ny282  - 8 I a. 

i=l 

N 

+ E log 
i=l  Y 


The  first  order  derivatives  are,  letting  y 


ia^G) /y 


(A2.5)  =»  log  L(6,G,T)  Ze  . z 

OP  P . _ 1 

1=1 


(A2.6)  5 log  L(8,G,Y)  = Ng  . 1 ” rh±l 

3G  y i=1  <J>(yi) 


(A2.7)  3 I°g.L(6’G'1f)  ■ -»e2V  - i “ y, 


3Y 


Y i=l 


Turning  to  the  second  order  derivatives  we  have 


( A2 . 8 ) 


32  log  L(8,G,y)  _ N _ 2 


36" 


= - F - NY/ 


(A2.9)  'j'~9B3G(B?G?y)  = N 


(A2.10)  92  '1afaL(8aG>y)  = "2ne^ 

. 383y 


<"•“>  s2-^1(l!-c'Tj ■-$  1 <*, * H^) 


(A2.13) 


92  log  L( 8 ,G i 


NB2  + -4  i X1 

y i=i 


(2-yK  ^Lr) 


Put  9 = ( B,G ,y ). Let  9 = (8  ,G  ,y  ) denote  the  likelihood  estimator 

and  0 = (0  ,G  ,y  ) the  true  parameter  set.  As  we  here  have  a regular 

o o o o 

likelihood  estimator  it  follows  that  the  limiting  distribution  of 
/N  (0-9q)  is  normal  with  expectation  zero  and  covariance  matrix 


(A2.ll*)  D2(^(0-0  ))  = (-  jp  E(-2-  "7s  /0=9  )) 

o N 90^  o 


where  3 ■ ^ given  by  eqs.  (A2.8)-(A2.13) . For  large  samples 

we  may  use  ^•9^/9=9)  computed  from  the  data  instead  of 

E(S  ^2^  L^9^/6=9°)  because 


Next,  consider  the  case  of  events  seen  jointly  by  two  stations. 
The  joint  distribution  of  loo  (A/T)  :s  a^  and  &2  i-s 


(A2.l6)  G(a^a^)  = f h^ ( a^/m)  exp( -dm)  dm 


/%dm2L£i>w^aife2i>)  „„(-8„)d. 


1/  + y'i 


2 2 
°2  + Y2 


1 a.-G. 


A(2.17)  h.  (a  ,/m)  = n — <I>(  — — 1)  exp[-(a.-B.-Q.-m)2/2a?] 

• 1 • C TT  O • Y # 1^3.  1 


u - a!  - BfS.  v=a2  - B2-  Qj 


a =*1  ‘ °1 


for  i = a,  2 


Then 


- 50  - 


(A2.l8 ) = / h;](a]/m/h0(a2/m)  exp(-Bm)  dir.  = 


- a r _L  exp[-(  (u-m)2/2a2+(v-m)2/2a2)]  exp(-Bn)  dm 

- y 2TOlaa 


Now 


(A2.19)  ( (u-m)2/ cr2  + (v-m)2/o2)-6m  = 


= 4 (ofu^v-ofofB))' 


12  1 2 


- 2(o2+gfT  ((u~v)2  + B ( 2a2u+2o2v-o2o2 B ) 


And  we  find 


(A2.20)  Fn  = 


Next 


$(_*)$(_£)  =4  exp[-((u-v)2+6’2a2u+2a2v-a2a2B))/2(a2+o2)] 

Yj  Y2  v^1t  U 2 


oo  m+^.+Q  -G.  m+B  +Q_-G^ 


(A2.21)  F2  = ; exp(-Bm) 

—oo  1 1 1 2 2 


dm 


Put  -VVV  C^2+G2  = d’sl  = ' S2  = ’/°FYI  ; 3 = /oi+a2+Yl+7I 


Then 


(A2.22)  F_  = / exp(-Rm)  dm  = 

^ oo  si  s2 


= -i-  exp(-Sm)]  + 

B S1  s2 


^ 


i 


- si  - 

cr> 

+ ~ / ♦(— — -)  JF-  exp[-(n-d)2/2s2]  exp f -Bm ) dm 

si  /2it  sp  2 

00 

1 , ^ ,m— 


— / $ (2^1  ysi  exp[-(m-c )2/2s? ] exn(-Bm) 
B s /2u  s,  1 


dm 


2 


(A2.23)  -(m-d)2/2s2  - 6m  = -(m-(d-s2B ) )2/2s2  — + ^ B2s2 


jimilarily 


(A2.2U)  -(m-c  )2 /2s2  - Bm  = -(m-(c-sJB)  )2/2s2  - eg  «■  B2s2 


So 


6“S  ” C c-s  *”d 

(A2.25)  = exp( -Bd+B2s2/2)  -■  $( — ) 1>  ^ ( ) exp(~Bc+B2s2/2) 

^ 2 6 /42+sz  6 /iTT-T  x 

S1  S2  * Sl+  2 


And  finally 


a-G,  G, 


(A2.26)  G(a  a ) = — )axo(-y/2) 

V2"VV°2  Yl  Y2 


(<t  ( 01 1 ) exp(a0)  +<S>(a,)  expfa. 


where 


(A2.2T)  y = ((a1-a2-(B1+Q1-B2-Q2))2+B(2a2(a1-B1-Qi)+2CT^(a2-B2-Q2)-a^CT2g))Aoc1+o:^ 


(A2.28)  a j = ((B,+Q1-B2-Q2)-(G1-G2)-(a2n2)B)y/pi+0247l+Y2 

(A2.29)  a2  = B(B2+Q2-G2)-*62(o2+Yp/2 

(A2.30)  a3  = ((B2+Q2-B1-Q1)-(G2'G 

(A2.31)  a4  = B(B1+Q1-G1)4B2(aJ+Y2)/2 


The  marginal  distrib’ition  of  a , given  that  the  event  is  seen  at  both 
stations  is 


4 


A2.32)  e(aL)  = / Gfa^J  da, , 


= ('o  / *{—-)  t-,:p[-((al-a^(B1+Qi-B,rQ2))2-28o2(ai ,-B, -Q2))/2(o^+o?J]  da? 


= C exp[-(B2a^-2B2o'4(a  -Bn-Qn  ))/2(a^+c2)]  ^ So’;+a^ 

O 1 lll'l  1^  If. 


VVarBrVBqj-G2;\ 


/'  O ■'■p 

. I / ' t- 
1>  I ( 


/q  L 2+Y? 


P'rom  this  w»-  can  obtain  the  eoridi  fionnJ  distribution  oi  a.,,  given  a j 
for  events  seen  at  both  stations.  We  get, 

(A2.33)  g(a2/aj)  = 0(a  2 )/g(a1  ) = 


* Ia-^)  ^ exp  [ - ( a2-a;L-B ) 2 /20  2 L 


where 


aj+B-G, 

* } 
2 


B * b2  ♦ VW60!  *»a 


n2  = n?  + o2 

1 2 


This  distribution  is  of  the  same  type  as  the  distribution  of  observed 

log  (A/T)  for  given  true  magnitude  eq.  (Al.l)  so  the  first  two  moments 

of  g(a.  /&  ) are  given  by  eq . (A1.10)  and  eq . (AI.l^).  On  the  other 
d 1 

hand,  had  w e also  considered  the  events  not  seen  at  station  2 given 
that  the  events  were  recorded  at  station  1 it  is  readily  seen  that  we 
btain  the  following  distribution 


(A2. 3h) 


g (a2/ai)daJ= 


°xp[  - ( a0-(  a +( -Q_  -o2R)  ) 2/?o2]  da 

Y2  V2II  a d L d 1 1 

' if  seen  at  station 


-(a,+(B2+Q2-B1~Q  -oJp))-G2 

<p  ( — - )if  not  seen  at  station  2 

,(  2 2 

r 
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The  log  likelihood  based  on  (A2.33)  for  a sample  of  N independent 
observations  is  except  for  a constant 


a2i-G 


(A2 . 35 )log  L=  Z-Nlogy  “(a2i*Q2i  “(ali_Qli)_B  } /2°  + log  ~ ) 

i ^2 

. ./aii+B+Q2i-Qli-°2, 

- ^ > 


where  B = (B^-B^-Bo2)  and  c2=o2+a 2 


Putting 

G = G^  and 

y = Y2  , V Q2i-Qu 

(A2.36) 

x.  = 
l 

(B+Q.+a1:i-G)/  l/a2+Y2 

(A2.3T) 

yi = 

(a2i-G)/Y 

and 

(A2.38) 

z(x) 

4>’  (x) 

$(x) 

The  first  order  derivatives  are 

N 


(A2.39) 


= z (a  _B-0  _a. . )/<J2-  , a z(x,  ) 


3B  . ' 2i  1 li ‘ 

1=1 


(A2.h0)  ^ *(yi)/Y  + ^7  *(*i) 

(A2.U)  = “f  (a2rB-Qi-aii)2/P^P;Y7yZ  (x.) 


(A2.U2)  -YflY^h  + (Wfr)  Z(xi} 


i-1 


The  likelihood  estimator  is  then  obtained  by  solving 

3 1 ogL  9 logL  3 logL  31ogL 

Tr  TG  ao  3 y 


0 


- $1.  - 

Lot  0 = ( li  ,o  ,y ) deriot'-  the  i.nn-  [>nriun'-l.cr  not.  rind  ft  = (|'./J,o,y) 

denote  the  likelihood  :;t  i mate . It.  then  follows  d i roof  I y from  t. tin- 
general  properties  of  the  maximum  likelihood  method  that  /n(0q-O) 
is  asymptotically  normally  distributed  with  expectation  zero  and 
covariance  matris  given  by 


(-  — E(d  w-°B-  ))-* 

1 N V 3e2/e=e0 ; ' 


which  may  be  estimated  by 

, I 32logL  1 
~ N 362/G=0 

Turning  to  these  second  order  derivatives  we  find  after  some 
calculation 

(Ar.1.3)  ^ 


(A2.UU) 


32logL 

3G2 


= E - 7.(y.)(y.+z(y.  ))/y2  + z(x  . ) (x+z(x.  ) ) 


i=l 


(A2.U5)  = E -(a2.-B-Qra2.)2.  3/a  + 

i=l 


+ f^7|y2  [a2(x,  (x,+z(x,  ))-2)  + Y2 ] 


1 1 


4 


i 

I 


(A  PM.) 

<)y 


N y,z(y.) 


>:  ■— t"  (r'-yi(yi+z(yi))  + 


+ [Y2(x.(xi+z(x.))-2)  + a2] 


(*2.1.7)  £ - ±4  <*i  * 

8BdG  . „ a^+vz 
i=l 


(A2'W)  ” - <*2i-B-V*,i)2/03- 

1 = 1 


oz(x. ) 

(aW)3/2  (xi(x.  + z(x.)-l) 


t N yz(x-.) 

(A2'h9)  = .Z=1  ~ (aV?) 3/2  (x.(xi+z(x.))-l) 


92logL 


N 


(A2.50)  (onY"T3/2  (^(x.+zCXi)-!) 


(A2l5l)  - V1- (yi(yi+z(yi)-l)  + 


+ ra^>)Y/?  (xi(xi+z(xi)  )-l) 


and  finally 


(A2.52) 

dO  o*Y 


N yox^zCx^^) 


ifi  - (02+V2)2  {3-xi(xi+2<*i))) 


